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AN EFFICIENT IMPLEMENTATION OF BREZZI-DOUGLAS-MARINI 
(RDM) MIXED FINITE ELEMENT METHOD IN MATLAB 

SHUN ZHANG 


Abstract. In this paper, a MATLAB package bdm_mf em for a linear Brezzi-Douglas- 
Marini (BDM) mixed finite element method is provided for the numerical solution of 
elliptic diffusion problems with mixed boundary conditions on unstructured grids. 
BDM basis functions defined by standard barycentric coordinates are used in the 
paper. Local and global edge ordering are treated carefully. MATLAB build-in 
functions and vectorizations are used to guarantee the erectness of the programs. 
The package is simple and efhcient, and can be easily adapted for more complicated 
edge-based finite element spaces. A numerical example is provided to illustrate the 
usage of the package. 


1. Introduction 

In recent years, MATLAB is widely used in the numerical simulation and is proved 
to be an excellent tool for academic educations. For example, Trefethen’s book on spec¬ 
tral methods m is extremely popular. In the area of finite element method, there are 
several papers on writing clear, short, and easily adapted MATLAB codes, for example 
[UEinnici]. Vectorizations are used in m and m to guarantee the effectiveness of the 
MATLAB finite elements codes. The mixed finite element [niEiii! is now widely used 
in many area of scientihc computation. For example, in [3 0 [3 019], we use RT(Rviart- 
Thomas)/BDM(Brezzi-Douglas-Marini) space to build recovery-based a posteriori error 
estimators. On the other side, except for the clear presentation of [2| on RTq, the imple¬ 
mentation of more complicated BDM elements is still somehow confusing for researchers 
and students. The purpose of this paper is to fill this gap by giving a simple, efficient, and 
easily adaptable MATLAB implementation of BDM\-Pq mixed finite element methods for 
of elliptic diffusion problems with mixed boundary conditions on unstructured grids. 

For linear BDM elements, there are several versions to write the basis functions ex¬ 
plicitly. In [3 [HI , basis functions are defined on each elements using heights and normal 
vectors. This version of basis functions is less straightforward than the basis functions 
defined by barycentric coordinates. After all, everyone is very familiar with barycentric 
coordinates in finite element programmings. Thus, in this implementation, we will use the 
definition which only uses barycentric coordinates. 

There are two basis functions on each edge for linear BDM elements. Unlike the RT 
element, BDM basis functions depend on the stating and terminal points of the edge. Thus, 
when assembling the local matrix, for a local edge on each element, we need to make sure we 
find its correct global stating and terminal points of the edge. There is an implementation 
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of BDM element in iFEM package [inj https://bitbucket.org/ifem/ifem/. But in 
order to make the local ordering of edges in an element is the same as the global ordering 
of edges, the triangles are not always counterclockwisely oriented. This will cause confusion 
for programmers. And if we use this ordering, sometime we may need two kinds of element 
map, one is counterclockwisely ordered, the other is ordered by the indices of vertices. 
This will make things more complicated. In this implementation, we will use the standard 
counterclockwisely ordering of vertices of triangles. 

Besides these issues, to get the right convergence order, we need to handle the boundary 
conditions carefully. In this package, we use the basic data structure of iFEM m, and 
full MATLAB vectorization is used. 

In a summery, in this package: 

(1) BDMi basis functions are explicitly defined using barycentric coordinates. 

(2) Elements are still positively oriented. A function is used to hnd the right global 
stating and terminal vertices of an edge in a local element. 

(3) Mixed type of boundary conditions are handled correctly to guarantee the right 
order of convergence. 

(4) MATLAB build-in functions and vectorizations are used to guarantee the erectness 
of the programs. 

The package can be download from http: //personal. cityu.edu.hk/~szhang26/bdm_ 
mfem.zip, 

The paper is organized as follows. Section 2 describes the model diffusion problem, 
its BDMi-Pq mixed finite element approximation and the corresponding matrix problem. 
Section 3 introduces a simple example problem to demonstrate out MATLAB code. In 
Section 4, edge-based linear BDM basis functions are defined. The main part of the code 
for the matrix problem is discussed in Section 6. MATLAB functions to checking the errors 
are discussed in Section 7. Finally, we discuss some related finite elements in Section 8. 


2. Model problem and BDMl-PO mixed finite element method 

2.1. Model problem. Let D be a bounded polygonal domain in 3?^, with boundary 
dVt. = F^ U Fjy, F^ n Fjy = 0, and measure (F d) / 0, and let n be the outward unit vector 
normal to the boundary. For a vector function r = (ti,T 2 ), define the divergence and 
curl operators by V • r = and V x r = For a function v, define the 

gradient and rotation operators by Vn = and V-*-n = (g^, ~g^)*- 

Consider diffusion equation 

(2.1) - V • (a(x)Vu) = / in D 
with boundary conditions 

(2.2) — aVtt • n = g'jv on Ftv and u = g£) on F^). 

We assume that the right-hand side / G L^(D), in and that g^ G L‘^{Tn), 

and that a(x) is a positive piecewise constant function. Define the flux by cr = —a{x)Vu, 
then we have 


(2.3) 


a = —Vu and V ■ a = f. 
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Define the standard H (div; D) spaces as 

H{div,n) := {t e :V-T e 

H]\f{div;Q) := {r G H{div; D) : r • n = 0 on Ttv}- 

Multiply the first equation in ( |2.3| ) by a r G HN{div; D) and integrating by parts, we get 

(a” V, r) = -(Vu, r) = (V • r, v) - {t ■ n, goho 

where (•, •)aj is the inner product on a domain uj. If w = D, we omit the subscript. Then 
the mixed variational formulation is to hnd (cr, u) G fl(div; D) x L^(D) with cr • n = gi\[ 
on Tjv, such that 


(2.4) 


(a ^cr, t) — (V • r, u) 
(V • cr, v) 


-{t ■\i,gD)vD V T G i?iv(div;D), 

(/, u) VuGL2(D). 


The existence, uniqueness, and stability results of (2.4) are well-known, and can be found 
in standard references, for example, [Ij. 


2.2. BDMi-Pq mixed formulation. Let T = {K} be a regular triangulation of the 
domain D. Denote the set of all edges of the triangulation by T := Sj U Sd ^ £n, where 
£i is the set of all interior element edges and £d and £m are the sets of all boundary 
edges belonging to the respective Td and Tat. For any element K gT, denote by Pk{K) 
the space of polynomials on K with total degree less than or equal to k. The iL(div; D) 
conforming Brezzi-Douglas-Marini (BDM) space [3j of the lowest order is defined by 

BDMi = {t :t\k G BDMi{K),yK gT} with BDMi{K) = Pi{Kf. 
and let BDMi^]\f = BDMi n (div; D). The piecewise constant space Pq is defined by 

Fh = {u : v\k G PoiK), VK G T}. 

For simplicity, we further assume that a is a positive constant in each element K G T- 
Let Ttv be the one dimensional mesh induced by 7” on F^r. Define 


Pi{Tn) = {v : v\e G Pi{E), 'iE G Tn]- 

Let gN,h be the projection of g^i on Pi (Tat)- The BDMi-Pq mixed finite element discrete 
problem then is to find {crh,Uh) G BDMi x Pq with cr^ • n = g^^h on Fat, such that 


(2.5) 


(a Va, ta) - (V • T/i, ua) = -(ta • n,g'A,)ro y th G BDMi^n, 

-(V • cta, Vh) = -if, Vh) y Vh G Pq. 


The discrete problem (2.5) has a unique solution, and the following error estimates hold 
assuming that the solution (cr, u) has enough regularity: 

(2.6) — a-h)\\o < Chf\\a~^^‘^a -\\2 and ||u — ua||o < ^(Fllulli-|-/i^||cr|| 2 ). 


where h is the mesh size and || • ||fc is the standard norm of Sobolev space H^. Thus, we 
should expect order 2 convergence of cr and order 1 convergence of u if the solution is 
smooth enough. 
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Let ctat G BDMi be the interpolation of on Fat (The construction of ctat will be 
explained in Section ?? in detail). Then cr^ = crj^ + ctq, with ctq g BDMi^n solves the 
following discrete problem: 

(2.7) 

J (a ^ao,Th) - {V ■ Th, Uh) = -{Th ■ n,gD)ro - ^o-p^,Th) ^ Th & BDMi^n, 

\-{V-ao,Vh) = -{f,Vh) + {V-aN,Vh) Mvh^Po- 


2.3. Matrix problem. Suppose that all edges are uniquely defined with fixed initial and 
terminal vertices, in Section [51 we will define two basis functions cji^ i and 0^ 2 for an edge 
T'j G T, j = 1, • • • , NE, with^i? is the number of edges. For simplicity, we assume that the 
total number of all edges in Sj and Sd is M, and BDMi^j^ = span{0j^^, 4'j,2d = 1 • • • , M} 
(The real mesh may has a different order of indices). The basis of Pq is very simple. For 
the element Kj, its basis is Ij which is 1 on Kj and 0 elsewhere. 

We want to compute the coefficient vectors x G of cr^ and y G IR^"*" of Uh with 

respect to the BDMi basis and Pq basis 

NE NT 

(2.8) cTh = '^ + XNE+j(t>j, 2 ) and Uh = '^ 

j=l k=l 


and 


(2.9) 


NE 

<Jn= ^ + XNE+j(t>j,2) 

j=M+l 


with Xj and x^E+j are determined by gN,h (discussed later in Section ??). Then 


M 


NT 




i=i 


k=l 


NE 


=-(ffD, • nn)ro - X] 0j,i) + a;j+ArB(a Vj,i,0j,2)); 

j=M+l 

M NT 

(xj(a“Vi,2: 0i,l) + Xj+NE{or^(t>i,2-,4>j,2)) “ y~!i/fc(0t.2^ l)Afc 

j=l k=l 

NE 

= -{go, </>i,2 • nn)ro - ^ {xj{a~^ 4 >i 2, 4 >j,i) + Xj+NEia~^(f^i,2^ ^j,2))\ 

j=M+l 

M 

~ X/ + Xj+NE{(f>j,2^ 

i=i 

NE 

= -{fA)Ke+ ^ {xjC^ ■ (t)j^i,l)Ke + Xj+neC^ ■ (f}j,2^'^)Ke) 
j=M+l 
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for i = 1, • • • , M, t = 1, • • • , NT. The coefficients Xj, Xj+ate, j = M + 1, - , NE are 

known from g^^h- Rewrite it as a matrix problem, we have 

( 2 . 10 ) 


C 0 


hi 

62 


where R is a 2M x 2M matrix and C is a NT x 2M matrix with entries: 

( 2 . 11 ) 

(^ ^ 2 , 1 ? ^j,l)? (^ ^ 2 , 1 ? 2)5 (q: ^ 2 , 2 ! 1)5 

Bi+M,j+M = («~^2,2) ^^^, 2 )) Q,J = -(V • l(), Cij+M = -(V • 1^)) 

where z = 1, • • • , M, j = 1, ■ ■ ■ , M, and £ = 1, • • • , NT; and the right hand side where 61 
is a 2M X 1 vectors and 62 is an NT x 1 vector with entries 

NE 

bli = -(5d, 02,1 • no)r£, - ^ Vi,!) 0j,i) + a^i+A^E(a“ Vi,!) 0j,2)); 

j=M+! 

NE 

bli+M = -(5D, 02,2 • no)ro “ X/ Vi,2) 0j,!) + 2^i+A^E(a“ Vi,2) 0j,2)); 

j=Af+! 

NE 

b 2 i = -{f,l)Ki+ ^ (xj,(V • 0 j 1, 1 )a:^ + Xj+ArE(V • 0 j 2, l)xj 

j=M+l 

where z = I,-- - ,M and ^ = I,-- - ,NT] (xi,-- - , xm, a^i+iVE, • • • , xm+ne, 2/!, • • • ^VntY 
is the 2M + NT x 1 solution vector. 


3. A Numerical Example 


We will demonstrate our MATLAB program by a simple test problem. Let 11 = (—1,1)^, 
with Ttv = {x G (—1,1)} X {y = 1}, and the rest is Tai The mesh is given in Fig. We 
choose the diffusion coefficient and the exact solution to be 


a = 


10 if X < 0; 
1 if X > 0; 


and m(x, y) 


(x^y^ + x)/10 + y if X < 0; 
x^y^ + X + y if X > 0. 


The right-hand side / = — 2(x^ + 2/^)- The exact cr is 

, \ / —(2xy^-|-1, 2x^y-|-10)* if x < 0; 

I —(2xy^-|-1, 2x^y-|-1)* if x > 0. 

It’s clear that a itself is not continuous, but its normal component is continuous. The 
boundary conditions are 


zz(x,y) 


' yVlO + y-1/10 
+ 2 / + 1 

(x^ -|- x)/10 — 1 
x^ -|- X — 1 

\ 


if X = — 1; 
if X = 1; 
if X < 0 and y 
if X > 0 and y 


- 1 ; 

- 1 ; 


and 


-aVu ■ (0,1)* = qn = 


—2x^ — 10 if X < 0 and y = 1; 
—2x^ — 1 if X > 0 and y = 1. 


We choose this exact solution to emphasize that the flux cr is in Lr(div;n), but not the 
Vzz. The main MATLAB codes of a, /, gn and yiv are given in exactalpha .m, f .m. 
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Listing 1. exactalpha, f, gD, and gN 


1 

function z 

= exactalpha(p) 

2 

ix = (p (: , 1)<0) ; z 

= ones(size (p,1),1); z (ix) = 10.0; 

3 

end 



4 

function z 

= f(p) 


5 

X = p ( : , 1) ; 

Y = P ( 

,2); z = -2*(x.*x+Y.*y); 

6 

end 



7 

function z 

= gD(p) 


8 

X = p ( : , 1) ; 

Y = P ( 

,2) ; 

9 

z= (x<0) . * (0 

.l*x.*x 

*Y■*Y+0.1*x+y)+(x>=0).*(x.*x.*y.*y+x+y); 

10 

end 



11 

function z 

= gN(p) 


12 

x = p ( : , 1) ; 

Y = P ( 

,2); z=(x<0) .*(-2 *x.*x-l0) + (x>=0) .*(-2*x.*x-l); 

13 

end 




gD . m, and gN. m respectively. Here the Dirichlet boundary condition is given by the exact 
solution in gD. m for simplicity. 



Figure 1. mesh of the example Figure 2. a triangle 

The MATLAB program to solve this problem is given in main.m. Lines 2-7 read the 
essential geometric data of the problem. Line 9 generates the edges, and other necessary 
geometric relations. Lines 11-12 solve the problem by the BDMi mixed finite element 
method. 



4. Triangulation and geometric data structures 

4.1. Geometric description and geometric relations. We follow m for the data 
representation of the set of all vertices, the regular triangulation T, the edges, and the 
boundaries. 

The set of all vertices AA = {zi, • • • , zat} is represented by an Nx2 matrix node (1 :N, 1:2), 
where N is the number of vertices, and the i-th row of node is the coordinates of the Lth 
vertex Zi = {xi,yi), node(i, :)= [xi,yi]. Lines 2-3 of main.m give the node of our 
example. For example, the 1st vertex has coordinates xi = —1 and yi = 1. 
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Listing 2. main.m 

1 %% Geometry setting 

2 node = [-1 1; 0 1; 1 1; -0.5 0.5; 0.5 0.5; -1 0; 0 0; 1 0; ... 

3 -0.5 -0.5; 0.5 -0.5;-1.0 -1.0; 0.0 -1.0;1.0 -1.0]; 

4 elem =[4 2 1;4 1 6;4 6 7;4 7 2;5 3 2;5 2 7;5 7 8;5 8 3;9 7 6;... 

5 9 6 11;9 11 12;9 12 7;10 8 7;10 7 12;10 12 13;10 13 8]; 

6 bdEdge = [2 0 0;1 0 0; 0 0 0;0 0 0;2 0 0;0 0 0;0 0 0;1 0 0;... 

7 0 0 0;1 0 0;1 0 0;0 0 0;0 0 0;0 0 0;1 0 0;1 0 0]; 

8 %% Geometric relations 

9 [edge,elem2edge,signedge] = geomrelations(elem); 

10 %% 

11 [sigma,u] = diffusionbdm(node,elem,bdEdge,elem2edge, edge, .. . 

12 signedge,@exactalpha,@f, SgD, @gN) ; 


Listing 3. geomrelations.m 


1 

function 

[edge,elem2edge,signedge] 

= geomrelations(elem) 


2 

NT = size(elem,1); 





3 

totalEdge 

= sort([elem(:,[2,3[ 

); elem(:,[3,1 

); elem(:, [1,2]) 

] , 2) ; 

4 

[edge, useless, j] = unique(totalEdge, 'rows 

); 


5 

elem2edge 

= reshape(j,NT,3); 





7 

signedge 

= ones(NT,3); 





8 

signedge( 

: , 1) = signedge(:,1) 

- 2* 

(elem(:,2 

>elem(:,3)); 


9 

signedge( 

:,2) = signedge (:,2) 

- 2* 

(elem(:,3 

>elem(:,!)); 


10 

signedge( 

:,3) = signedge (:,3) 

- 2* 

(elem(:,1 

>elem(:,2)) ; 


11 

end 







The triangulation T is represented by an NT x 3 matrix elem(l:NT,l:3) with NT the 
number of elements. The i-th element Ki = convjzj, Zj, z^} is stored as elem (i, : ) = ... 
[i j k], where the vertices are given in the counterclockwise order. Lines 4-5 of main.m 
give the elem of our example. For example, the 1st element Ki has three vertices in the 
order of Z 4 , Z 2 and zi. 

We call the the opposite edge of the i-th vertex, i = 1, 2, 3 of a triangle the z-th edge of 
the triangle. 

The matrix bdEdge (1: NT, 1:3) indicates which edge of an element is on the boundary 
of the domain. For a non-boundary edge, the value is 0; the value is 1 or 2 for a Dirichlet 
or Neumann boundary edge respectively. Lines 4-5 of main.m give the bdEdge of our 
example. For example, the 1st element Ki has three edge, the first edge is on the Neumann 
boundary, and the other two are not on the boundary. 

The MATLAB code geomrelations .m generates edge, elem2edge, and signedge. 
The explanations of these codes can be found in jini . 

The matrix edge (1 : NE, 1 : 2 ) is dehned with NE the total number of edges, and the A:-th 
edge Ek with the starting vertex zi and the terminal vertex Zj is stored as edge (k, :) = ... 
[1 j]. We always ensure that edge (k, 1 ) < edge (k, 2 ). Lines 3-4 of geomrelations .m 
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Listing 4. gradlambda.m 



generate edge. For our example, NE=2 8 and edge= [1 2; 1 4; 1 6; 2 3; 2 4; 2 ... 
5; 27; 35; 38; 4 6; 47; 57; 58; 67; 69; 611; 78; 7 9; 710; .. 

7 12; 8 10; 8 13; 9 11; 9 12; 10 12; 10 13; 11 12; 12 13] . For example, the 

2 nd edge’s starting vertex is zi and terminal vertex is Z 4 . 

For an edge Ej with the starting vertex Zg and the terminal vertex zt, its normal 
direction is defined as {yt — ys,Xs — xtY/\Ej\. 

The matrix elem2edge (1 :NT, 1 : 3) is the matrix whose A:-th row represents the 3 
global indices of the edges in the order of local edges. Line 5 of geomrelations .m gen¬ 
erates el em2 edge. For our example, elem2edge= [12 5; 3 10 2; 14 11 10; 7 ... 

5 11; 4 6 8; 7 12 6; 17 13 12; 9 8 13; 14 15 18; 16 23 15; 27 24 23; ... 

20 18 24; 17 19 21; 20 25 19; 28 26 25; 22 21 26]. For example, the 2nd ele¬ 

ment’s three edges are E^, Eiq, and E 2 . 

Each edge has its fixed global orientation, while on each element, it has a local orienta¬ 
tion, localEdge= [2 3; 3 1;1 2]. The matrix signedge is an NT X 3 matrix with value 1 
denoting the local and global orientations are the same, and —1 denoting that they are dif¬ 
ferent. For our example, signedge=[-l 1 - 1 ; 1 -1 - 1 ; 1 -1 1; -1 1 1 ; -1 ... 

1 - 1 ; 1 -1 - 1 ; 1 -1 1 ; -1 1 1 ; -1 1 - 1 ; 1 -1 - 1 ; 1 -1 1 ; -1 11 ; -11 .. 
- 1 ; 1 -1 - 1 ; 1 -1 1 ; -1 1 1 ]. 


4.2. Barycentric coordinate and its gradient. On an element K with counterclock¬ 
wise vertices {zi, Z 2 , Z 3 }, we define the barycentric coordinate A; = {uiX + biy + Ci)/{2\K\), 
i = 1,2,3, such that Xi{zj) = 5ij. We only need to compute the gradient (and curl, rot, 
div) operators in this paper, so we only need to compute a;, bi, and the area \K\. The 
formulas are 


VA; 


1 

Wl 


Ui+i - yi+2 \ 

7 ^ 1+2 J 


and 2 |iL| = det 


X2 - Xl X3 - Xl 
y2 - 2/1 2/3 - 2/1 


The function gradlambda.m computes the coefficients a, b, and area. Here a and b are 
two NTx3 matrices, with each row stores the coefficients a and b for the three barycentric 
coordinates of corresponding 3 vertices. 


4.3. Normal and tangential vectors. On a local element K with counterclockwise 
oriented vertices {zi, Z 2 , Z 3 }, and Ei = {z 2 ,Z 3 }, E 2 = {z 3 ,zi}, and E 3 = {zi,Z 2 } (Figure 
[^, we will discuss the V and V"*- of A 2 as examples: 


V^Aa 



and V'''A 2 



2/3 - 2/1 
Xl - X3 


Xl - X3 
2/1 - 2/3 
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where \K\ is the area of K. On edge Ei = {z 2 -, 2 : 3 }, the unit tangential and normal vectors 
are 


tsi = 


1 


|£^il 


^3 - X2 

ys - 2/2 


and n^i = 


1 


where |-E| is the length of E. From Figure 


we have 


l^^il 

Xl 

yi 


2/3 - 2/2 

X2 - X3 
X3 
2/3 


• n^i = 


2/3 

Xl 


2/1 

3^3 


tEi = —Hi, with Hi is the height of Ki onEi. Now, by the fact 2\K\ = |Fli| • Hi, we have 

V-*-A2 • = VA 2 • = —l/|Fli|. 

Similarly, 

V-*“A3 • UEi = VA3 • = 1 /l-Eil- 

On the hand, since the tangential and normal vectors of an edge are orthogonal, we have 
(4.1) V-^A2 • ngj = V-^A3 • = 0. 

Globally, on E^ = {zs,zt} {s <t), the unit tangential and normal vectors are 


tEo = 


1 


Ep 


Xt - Xs 

yt - ys 


and ng, = 


1 


Ep 


yt - ys 

Xs - Xt 


We call the adjacent element whose out unit normal vector is the same as as , and 
the element whose out unit normal vector is the opposite of as K^. On both and 
, we have 

(4.2) V-^As • = VAs • =-I/IFI^I and V-^A* • = VA* • = 1/I.E^I. 


5. Constructions of edge-based BDM basis functions 


Let Aj be the standard linear Lagrange finite element basis function of the vertex Zi, i.e., 
it is piecewise linear on each element, globally continuous, 1 at Zi and 0 at other vertices. 

For an edge Ep = {zg, zt}, s < t, 1 < i < NE with the globally fixed starting vertex s 
and terminal vertex t, its two BDMi basis functions associated with the edge are 

(5.1) cjip^i = XgV^Xt and (f)p ^2 = —AtV^’^A^. 

It’s clear that the basis functions are only non-zero in the two adjacent elements and 
(one element in the case of boundary) of Ep. 


Lemma 5.1. There hold 

f 0 if k £■, 

(5.2) (f>p^i-nEi^\Ep. — I Xs/\Ep\ if k = £' ^^,2'|E*. 


0, if k ^ £■, 

Xt/\Ep\, ifk = £. 


and 


(5.3) V • cfp p = V • (j)p^2 = 2\K-\ ^ V • cj)p^^ = V • = ~ 2\K+\ 

Proof. By the discussion in Section [4.3t we have 

X/-^Xs ■ UEi =-l/\Ep\ and Xf nEi = l/\Ep\, 


on . 


and 


4>p i ■ he^Iei = ^s/\Ep\ and 4>p^2 ' ^Ei\ep = \l\Ep 


Assume ’s three counterclockwisely ordered vertices are {zr-, Zg, zt}■ We call edge 
with two endpoints zt and Zr- to be Ep^r-, and the edge with two endpoint Zr- and Zg 
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to be Eg^r-- Here the orientations of these two edges are not important. Then by (4.1), 
V-*-As • = 0, and the fact A* is zero on Eg^r-, we get 

cf)f^ 2 ■ = 0, with E = Eg^r- or Et^r-- 

We can prove (5.2) is true for (pn and similarly. 


The divergence part of the theorem can be easily proved by recalling the definition of 
Xg and Xt on and K^, and notice that the counterclockwise ordering of vertices on 
is {zr-, Zg, zt}, and that on is {zr+, zt, Zg}. □ 

Remark 5.2. There are other definitions of the BDMi basis functions. One of the popular 
choice is the following: 

(5.4) = AsV-^At - AiV-^As and (/>£ 2 = XgV^Xt + XtV^Xg. 

The first basis function coincides with the RTq basis function with property (fn ■ ue^. = 


\Ei\5ik. This set of choices of basis functions is hierarchical. The reason we choose (5.1) 
is that it is symmetric. If the reader needs the basis to be hierarchical, one should choose 


(5.4) 


Remark 5.3. The other versions of basis functions, for example, the RT basis function 
used in [2], choose (jTf*' • ngj, = 5ke- The issue of this choice of basis functions is that the 
mass matrix has a condition number C/imax/^min? where hmax cind hmin o,re the maximal 
and minimal diameters of the elements respectably. This will be a problem for an adaptively 
generated meth. Though it can be easily fixed by preconditioning it with the inverse of the 


diagonal matrix of it, we avoid it by choosing basis in (5.1) or (5.4). 


6. Constructing and solving the matrix problem 


6.1. Assembling matrices. 

6.1.1. Local matrices. For an element K, we want to compute the local contributions of 
this element to the matrices B and C. 

We use localEdge= [2 3; 3 1; 1 2 ] to denote the local edge with respect to the local 
indices of vertices. For the i-th edge, we let ill = localEdge (i, 1) and 112 = ... 
localEdge (1,2) to denote the local starting and terminal vertices of it. When the local 
orientation and the global orientation of the edge E is different (slgnedge of the edge is 
— 1), we need to switch the local order to find the right global starting and terminal vertices 
of edge. For a given i-th local edge, by the line 11 = (slgnedge (:, 1) >0) . *111+ . . . 
(slgnedge (:, 1) <0) . *112, we get 11=111 if the local orientation and the global orien¬ 
tation are the same, and 11=112 otherwise. i2 can be done similarly. Once we find the il 
(local starting vertex) and i2 (local terminal vertex) with right global orientation, we can 
get the corresponding coefficients an and hn of Aji = {anx + buy + Cji)/(2|A|), and the 
same things for i2. The function BDMrlghtorder .m does the above job. With an input of 
i = 1, 2, or 3 be the local vertex index of an element, this function returns the values 11 
and 12, which are the correct starting and terminal vertices of the corresponding edge with 
respect to the global fixed edge orientation, and all,al2,bll,bl2 are the corresponding 
coefficients of Aji and Aj 2 - 

To compute the local contribution of an element K to B, the local mass matrix contains 
three cases, (a“ Vi,i, ^j,i), («“ Vi,i, and (a“ with 

^i,i = AjiV'''Aj2, 4 >i ,2 = —•^i2V'''Aji, = XjiV~^Xj2, and 0^2 = —Xj2V~^Xj2. 
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Listing 5. BDMrightorder.m 



For barycentric coordinates, fj^XiXjdx = (1 + 6ij)\K\/12. An easy computation shows 
that 

= a~^{l + 6iiji){ai2-aj2 + bi2-bj2)/i‘i8\K\)-, 

+ (^iij2)(ai2 • Oji + ^*2 • ^ji)/(48|Ar|); 

{a~^ct}i^2:4>j,2)K = a~^{l + 6i2,j2)iaii-aji + bii-bji)/{4:8\K\). 

For A; = 1, 2, V • (/)j ^ = 1/(2|A'|) when the Lth edge has the same orientation as the global 
edge, and —l/(2|iL|) otherwise. Thus, denote s{i) be the value of signedge (: , i) , we 
have 

-(V • 1)a: = -(V • = -s(*)/2- 

The local matrix (here we abuse the notations to use local * = 1---3, k = 1,2 to 
denote the local BDM basis functions): 

— ((V0;^^1, 1)k, (V02,1) 1)a'; (V03 1, 1)a'; (^<^3 2) K: (V<^2,2) ^)k-, (V03 2j 1)a') 

is simply —(s(l), s(2), s(3), s(l), s(2), s(3))/2, or - [signedge (:, 1: 3) , signedge { :, 1 : 3) ] 
in MATLAB code. 

6.1.2. Assembling global matrices. For a local index i, its corresponding global edge index 
is double {elein2edge {:, i) ) • MATLAB function sparse is used to generate the global 
matrices from local contributions. This is one of the key step to ensure the vectorization 
of the MATLAB finite element code, see mm for more detailed discussions on sparse. 
Here are some comments of the MATLAB code assemblebdm.m. 

• Line 1: The function is called by 

A = assemblebdm(NT,NE,a,b,area,elem2edge,signedge, inva) 
where A is a (2*NE+NT) * (2*NE+NT) matrix, a, b, area are the coefficients of 
local barycentric coordinates, signedge is the NT*3 matrix about the local and 
global edge orientations, and inva is NT*1 vector of a~^. 

• Lines 4-15: We generates the matrix B = {a~^crh,Th), (Th and Th in BDMi by 
assembling local element-wise contributions. 

• Lines 6-7: We generates the right starting and terminal indices of a global edge in 
the local element, and their corresponding coefficients of local barycentric coordi¬ 
nates. 

• Lines 8-10: We compute E, H, and G, which are <^^•^ 2 )^', 

and 2 ) ^j, 2 )A') respectively. 

• Lines 11-13: B is generated by sparse. 
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Listing 6. assemblebdm.m 


1 function A = assemblebdm(NT,NE,a,b,area,elem2edge,signedge,inva) 

2 B = sparse (2*NE, 2*NE); C = sparse (NT, 2*NE); D = sparse(NT,NT); 

3 %% Blocal (6,6) = [E(3,3) H(3,3); H'(3,3) G(3,3)] 

4 for i = 1:3 

5 for j = 1:3 

6 [il,i2,ail,bil,ai2,bi2] = BDMrightorder(i,signedge,NT,a,b); 

7 [j1,j2,aj1,bj1,aj2,bj2] = BDMrightorder(j,signedge,NT,a,b); 

8 E = inva./(48*area) .*(1+ (il==j1)) .* (ai2.*aj2+bi2.*bj2); 

9 H = -inva./(48*area).*(1+(il==j2)).*(ai2.*ajl+bi2.*bj1); 

G = inva./(48*area).*(1+(i2==j2)).*(ail.*aj1+bil.*bj1); 
ii = double(elem2edge(:,i)); jj = double(elem2edge(:,j)); 

B = B + sparse (ii, j j , E, 2 *ne, 2jcNE) Tsparse (ii, j j+NE, H, 2*NE, 2*NE) .. . 

+ sparse ( j j+NE,ii,H,2*NE,2*NE)+sparse(NE + ii,NE+jj,G,2*NE,2*NE); 


1:3 


C + sparse(1:NT,ii,-signedge(:,i)/2,NT,2 *NE) . .. 
+ sparse(l:NT, ii+NE,-signedge(:,i)/2,NT,2*NE); 

[B C';C D]; 


10 



11 



12 



13 



14 


end 

15 

end 


16 

for 

i = 

17 


ii ^ 

18 


C = 

19 



20 

end 


21 

A = 

[B 1 

22 

end 



• Lines 16-20: C is generated by local contributions [-signedge (:, 1:3) , ... 
-signedge(:,1:3)] . 

• Lines 21: A is generated by adding a zero matrix D on 22 block. 

6.2. Assembling the force / term. Since V ■ cr^ € Pq, we only need the numerical 
integration of the / term to be accurate as if / is a constant on each element. Thus, the 
term related to — (/, 1)k can be computed by a one-point quadrature rule: 

/ /dx f{^Xinidiymid)\P\i 

Jk 

where {xmid, ymid) are the coordinates of the gravity center of the element K. 

6.3. Generating boundary data. On the boundary, we need sign_D and sign_N to 
denote the difference between orientations inherited from the local edge ordering of the 
element and the global edge orientations like edge sign. Since the unit out normal vector 
of an element on the boundary is the same as the unit out normal vector of the whole 
domain, sign_D and sign_N are actually the difference of normal directions of Dirichlet 
and Neumann edges and global out normal directions. 

The MATLAB function boundary.m generates the Dirichlet and Neumann edges and 
their orientations with respect to the global edges. 

• Lines 3-10: We generates un-sorted Dirichlet and Neumann edges and their signs. 

• Lines 11-13: We generates sorted Dirichlet and Neumann edges and their signs. 

We denote the set of indices of edges on the Dirichlet and Neumann boundary to be 
ind_D and ind_N, respectively. 
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Listing 7. boundary, m 

1 function [Dirichlet, Neumann, sign_D, sign_N] = boundary(elem,bdEdge) 

2 NT = size(elem,1); 

3 totalEdge = [elem(:, [2,3]); elem(:, [3, 1 ]); elem(:, [ 1,2]) ] ; 

4 IsBdEdge = reshape(bdEdge,3*NT,1); 

5 Dirichlet = totalEdge((IsBdEdge == 1),:); 

6 Neumann = totalEdge((IsBdEdge == 2),:); 

7 NE_D = size (Dirichlet,1); NE_N = size(Neumann,1); 

8 sign_D = ones(NE_D,1); sign_N = ones(NE_N,1); 

9 sign_D = sign_D(:,l) - 2* (Dirichlet(:,1)>Dirichlet (:,2)); 

10 sign_N = sign_N(:,l) - 2* (Neumann(:,1)>Neumann(:,2)); 

11 Dirichlet = sort(Dirichlet,2); Neumann = sort(Neumann,2); 

12 [Dirichlet, I_d] = sortrows(Dirichlet); [Neumann, I_n] = sortrows(Neumann); 

13 sign_D = sign_D (I_d) ; sign_N = sign_N (I_n) ; 

14 end 


To computer the term — (r • n, 5 (p))r^, we need to be careful about two things. One 
is the difference between unit out normal vector of the domain and that of the edge, the 
other one is the numerical quadrature formula. In order to guarantee the convergence 
order of BDMi element, we should use two-point numerical quadrature on an edge. The 
2-point Gauss-Legendre quadrature of a function f{x) on interval [a, b] is: 


( 6 . 1 ) 


f{x)dx 


b — a 


a + b b — a 

2 V 3 


+ / 


a + b b — a 
2 ^ 2\/3 


Note that on a Dirichlet edge Ej = {zs,zt}, s < t, cf)j i ■ n£u = \s/\Ej\ and cj>j 2 • = 

Xt/\Ej\. Thus, to compute — X4>j,e''’^n) 9 Ddx by formula (6.1), we need the value of qd 

ni +n2 n2 -ni 


at quadrature points pi = 


and p 2 = 


nin2 ^ n2 - ni 


where ni and 


2 2y/3 2 2 V 3 

n 2 are the coordinates of the Zg and zt, respectively. We also need to know the value of 
\s and Xt at pi and p 2 , which are 


A.(pi) = Mn) = ^ ^ 


and Xsip 2 ) = Xtipi) = ^ - 


Thus, by letting s{j) = sign_D ( j) , s(j) = 1 or —1 be the number denoting the difference 
between the local and global edge orientation on Tp), we have 


- (^j,i • nn)gDdx 


- / • T^u)gDdx 

J Eq 


-s{j) (»(Pi)(l + ^) + 9D(n)(\ - ^)) , 

-s{i) (»(pi)(l - + + j^)) ■ 


To handle the Neumann boundary condition, on each Ej = {zs,zt} € £n, we want to 
compute Then cr^y • = Cj^i(f)j i ■ -j- Cj^ 2 (t>j ,2 ' Tet 

s{j) = 1 or —1 be the number denoting the difference between the local and global edge 
orientation on Ej, we should have 


s{j){cj,iXs + Cj^2h)/\Ej\ Ri qn- 
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Listing 8. rhside.m 

1 function [A,b,sol,freeDof]=rhside(node,elem,edge,bdEdge,area,A,sol,f,gD,gN) 

2 NT = size(elem,1); NE = size(edge,1); 

3 [Dirichlet, Neumann,sign_D, sign_N] = boundary(elem,bdEdge); 

4 %% Assemble right hand side. 

5 mid = (node(elem(:,1),:)+node(elem(:,2),:)Inode(elem(:,3),:))/3; 

6 b2 = accumarray((1:NT)',-f(mid).*area,[NT 1]); 

7 %% Drichelet BC 

8 nl= node(Dirichlet(:,1),:); n2= node(Dirichlet(:,2),:); 

9 pl=(nl-n2)/2 *sqrt(l/3) + (n2+nl)/2; p2=(n2-nl)/2*sqrt(l/3) + (n2+nl)/2; 

10 intgDphiln = (gD(pi)*(1 + sqrt(1/3))+gD(p2)*(1-sqrt(1/3) ))/4 ; 

11 intgDphi2n = (gD(p2)*(1 + sqrt(1/3))+gD(pi)*(1-sqrt(1/3)))/4 ; 

12 [useless, ind_D] = intersect(edge, Dirichlet, 'rows'); 

13 bbl = accumarray(ind_D,-sign_D.*intgDphiln, [NE 1]); 

14 bb2 = accumarray(ind_D,-sign_D.*intgDphi2n, [NE 1]); 

15 bl = [bbl;bb2]; 

16 %% Neumann BC 

17 nl= node(Neumann(:,1),:); n2= node(Neumann(:,2),;); 

18 pl=(nl-n2)/2 *sqrt(l/3) + (n2+nl)/2; p2=(n2-nl)/2*sqrt(1/3) + (n2+nl)/2; 

19 edgeLength_N = sqrt(sum( (nl-n2) .'2,2) ); 

20 intgNlams = edgeLength_N.*(gN(pi)*(1+sqrt(1/3))+gN(p2)*(1-sqrt(1/3)))/4; 

21 intgNlamt = edgeLength_N.*(gN(p2)*(1 + sqrt(1/3))+gN(pi)*(1-sqrt(1/3)))/4 ; 

22 [useless, ind_N] = intersect(edge, Neumann, 'rows'); 

23 sol(ind_N) = 2*sign_N.*(2*intgNlams-intgNlamt); 

24 sol(ind_N+NE) = 2*sign_N.*(2*intgNlamt-intgNlams) ; 

25 %% modify right hand side 

26 b = [bl;b2]; b = b - A*sol; 

27 %%freeDof 

28 isBdDof = false (2*NE+NT,1); isBdDof(ind_N)=true; isBdDof(ind_N+NE)=true; 

29 freeDof = find(-isBdDof); 

30 end 


Let the L^-projection of on to Pi{Ej) = spanjA^, A^} to be gi:<[,h = dgXs + dtXt- 
V {>^sAt)Ej (AtjAt)^;^ J \ J V i9NAt)Ej )' 


Since J^XiXjds = \E\{1 + 6ij)/6, replace {gN,Xs)Ej and {gN,Xt)Ej, by Ij^s and Ij^t using 
the two-point numerical quadrature (6.1) as we did for g£), we get 

ds = - l 2 j,t)/\Ej\ and dt = - ‘^Ij,s)/\Ej\. 

So ctnIej = + Cj, 2 (l>j ,2 with Cj^i = s{j){4:Ij^s - 2Ij^t) and Cj ,2 = - 2/^,^), 

or, in the notation of (2.9), 

Xj = s(j)(44s - 2Ij^t) and xne+j = s(j)(4/j,t - 2/^,^). 


With this know a tv , then the right-hand side of the discrete problem can be easily handled 
as in [U [To] . 

The followings are some comments about the function rhside .m of handling the / term 
and boundary conditions. 

• Lines 5-6; We generate terms related to — (/, u). 
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Listing 9. diffusionbdm.m 

1 function [sigma,u] = diffusionbdm(node,elem,bdEdge, elem2edge, edge, .. . 

2 signedge,exactalpha,f,gD,gN) 

3 NT = size(elem,1); NE = size(edge, 1) ; 

4 sol = zeros(2*NE+NT,1) ; 

5 inva =1./exactalpha((node(elem(:,1))+node(elem(:,2))+node(elem(:,3)))/3); 

6 [a,b,area] = gradlambda(node,elem); 

7 A = assemblebdm(NT,NE,a,b,area,elem2edge,signedge,inva); 

8 [A,b,sol,freeDof] = rhside(node, elem,edge,bdEdge,area,A,sol,f,gD,gN); 

9 sol(freeDof) = A(freeDof,freeDof)\b(freeDof); 

10 sigma = sol(1:2*NE); u = sol(2*NE+1:end); 

11 end 


• Lines 8-15: We generate terms related to — (r • n,g'£))r^. 

• Lines 17-24: We construct the ctat. 

• Line 26: We handle the right-hand side of the matrix problem. 

• Lines 28-29: We find the degrees of freedom of the matrix problem. (Those terms 
related to a-jy are known, thus not free.) 

6.4. Solving the BDM mixed problem. The MATLAB function diffusiondbm.m is 
the principle function to solve the BDM mixed problem. With all the building blocks 
given before, it is relatively easy. Line 5 is used to got a~^ in each element. All the rest 
lines are self-explanatory. 


7. Checking errors 

The hnal step of a numerical test is often the convergence test by computing some 
norms of the error between the exact and numerical solutions obtained for different mesh 
sizes. In our case, we need to compute 

||a“^/^(cr - cr/i)||o and \\u - Uh\\o. 

For the both terms — (t/i)||o and ||m —u/i||o, the direct way to compute them is 

summing up the errors on each element by direct computations. To lower the numerical 
quadrature order, the first step should be 

||Q;“^/^(cr - (Th)\\l = {a~^{cr - CTh),(T - a^) = {a~^a,a) - 2{a~^a,ah) + (a“V/j,cr,,), 

\W-Uh\\l = {u-Uh,u-Uh) = {u,u)-2{u,Uh) + {uh,Uh). 

The terms cr) and (u, u) can be computed exactly by softwares like Mathematica 

and Maple. In our example, {a~^cr,cr) = 1993/75 « 26.5733 and {u,u) = 18131/7500 R3 
2.41747. The terms {a~^crh, cTh) and {uh, Uh) can also be computed exactly. For the terms 
{a~^a,ah) and {uh,Uh), numerical quadratures must be used. 

Since this part of codes is less important, we only give brief comments about the func¬ 
tions we used. The MATLAB function exact s igma .m returns the exact value of cr. The 
function sigmahBDM.m returns the values of Ch at those numerical quadrature points. The 
function errorBDM.m is used to compute the errors. Here, we use a 6-point quadrature to 
compute them. Matrix xw is a 6 * 3 matrix, with xw (:, 1: 2 ) are the coordinates such that 
the quadrature point is p=pl* (1-xw (i, 1) -xw (i, 2) ) +p2*xw (i, 1) +p3*xw (i, 2) , where 
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Listing 10. exactsigma.m 

1 

function [sigmal,sigma2 

] = exactsigraa(p) 

2 

X = p ( : , 1) ; y = p ( : , 2 ) ; 


3 

sigmal=-2 *x.* y.* y-l; 


4 

sigma2=(x<0).*(-2*x.*x. 

*y-10 ) + (x>=0) . (-2*x.*x.Jcy-l) ; 

5 

end 



Listing 11. sigmahBDM.m 

1 function [sigmahl,sigmah2] = ... 

2 sigmahBDM(elem,node,NE,eleraZedge,signedge,wl,w2,sigma) 

3 [a,b,area] = gradlambda(node,elem); NT = size(elem,1); 

4 lambda_at_w = [l-wl-w2,wl,w2]; 

5 sigmahl = zeros (NT,1); sigmah2 = zeros(NT,1); 

6 for i = 1:3 

7 [il,i2,ail,bil,ai2,bi2] = BDMrightorder(i,signedge,NT,a,b); 

8 ii = double(elem2edge(:,i)); 

9 sigmahl = sigmahl+sigma(ii).*lambda_at_w(il)'.*bi2./area/2 ... 

10 -sigma(ii+NE) .*lambda_at_w (i2) ' .*bil./area/2,• 

11 sigmah2 = sigmah2-sigma(ii).*lambda_at_w(il)'.*ai2./area/2 ... 

12 +sigma(ii+NE) .*lambda_at_w (i2) ' . *ail./area/2,• 

13 end 

14 end 


pi, p2, and p3 are the coordinates of three vertices; xw (:, 3) is the corresponding weight 
for the point. 

We use the code listed in Listing to compute the error. 

The original mesh given has mesh size h = 1. We refine the mesh uniformly several 
times, for example, using [node, elem, bdEdge] =uniformbisect (node, elem, bdEdge) , 
where uniformbisect is a uniform refinement MATLAB function given in the package 
zFEM |10] . We compute the errors for different mesh sizes, and have the Table We get 

||a~^(<T - crfe)||o _ \\u-Uh\\o _ 

\\a-^{a - a-h/2)h'^ \\u - Uh/2\\o ^ 


This result is in agreement with the a priori error esteems (2.6). 


8. Related finite elements 

Once we know how to programming BDMi methods, we can actually write codes for 
similar finite elements, with two things in mind: Writing basis functions in barycentric 
coordinates and correcting the local edge orientation. For more higher order elements, we 
may use Legendre polynomials instead of Lagrange polynomials, but the idea is similar. 

8.1. Other H(div)-conforming edge-based basis in 2D. For RTo element, as men¬ 
tioned in Remark |5.2[ its basis function on an edge can be written as 

- AtV^A, 
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Listing 12. errorBDM.m 

1 function [error_sigma, error_u] = errorBDM(node,elem,NE,area,elera2edge, ... 

2 signedge,inva,sigma,u,ss,uu) 

3 nl = elem(:,l); n2 = elem(:,2); n3 = elem(:,3); 

4 pi = node(nl,:); p2 = node(n2,:); p3 = node(n3,:); 

5 xw=[0.44594849091597 0.44594849091597 0.22338158967801;... 

6 0.44594849091597 0.10810301816807 0.22338158967801;... 

7 0.10810301816807 0.44594849091597 0.22338158967801;... 

8 0.09157621350977 0.09157621350977 0.10995174365532;... 

9 0.09157621350977 0.81684757298046 0.10995174365532;... 

10 0.81684757298046 0.09157621350977 0.10995174365532]; 

11 NP =size(xw,l); NT = size(elem,1); 

12 uuhlocal = zeros (NT,1); % (u, u_h)_K 

13 sshlocal = zeros(NT, 1); % (\a'" {-l}sigma, sigma_h) _K 

14 shshlocal = zeros(NT,1); % (\a"" {-l}sigma_h, sigma_h) _K 

15 for i=l:NP 

16 p=pl* (1-xw (i, 1) -xw (i, 2 ) ) +p2 *xw (i, 1) +p3^cxw (i, 2 ) ; 

17 uuhlocal = uuhlocal + exactu(p)*xw(1,3); 

18 [sigmal,sigma2] = exactsigraa(p); 

19 [sigmahl,sigmah2] = ... 

20 sigmahBDM(elem,node,NE,elem2edge,signedge,xw(1,1),xw(1,2),sigma); 

21 sshlocal = sshlocal + (sigmal.*sigmahl+sigma2.*sigmah2)*xw(1,3); 

22 shshlocal = shshlocal + (sigmahl.*sigmahl + sigraah2.*sigmah2)*xw(1, 3) ; 

23 end 

24 uuhlocal = u.*uuhlocal.*area; uhuhlocal = u.*u.*area; 

25 sshlocal = inva.*sshlocal.*area; shshlocal = inva.*shshlocal.*area; 

26 error_u = sqrt(abs(uu-2*sum(uuhlocal)+sum(uhuhlocal))); 

27 error_sigma = sqrt(abs(ss-2*sum(sshlocal)+sum(shshlocal))); 

28 end 


Listing 13. checking errors 

1 uu=18131/7500; ss=1993/75; 

2 [error_sigma, error_u] = errorBDM(node,elem,NE,area,elem2edge, ... 

3 signedge, inva,sigma,u,ss,uu); 


Table 1. Errors of a and u on different mesh sizes 


h 

O 

b 

1 

7 

\\a ^(cr - cr;,)||o 

k “ 0 

\\u - UhWo 

||a-i(cr - (Th/2)\\o 

\\u - Uh/2\\o 

1 

1.6968e-01 


4.9712e-01 


1/2 

4.2091e-02 

4.0314 

2.4400e-01 

2.0374 

1/4 

1.0600e-02 

3.9707 

1.2118e-01 

2.0135 

1/8 

2.6630e-03 

3.9805 

6.0481e-02 

2.0037 

1/16 

6.6739e-04 

3.9901 

3.0226e-03 

2.0009 

1/32 

1.6705e-04 

3.9952 

1.5111e-03 

2.0002 

1/64 

4.1788e-05 

3.9975 

7.5555e-04 

2.0001 
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Actually, it is insensitive to the staring and terminal vertices: 


cfif = A.V^At - AtV^A, = AtV^A, - A,V^At, 

and on adjacent elements K~ = {zr-, Zg, zt} and = {zr+, zt, Zg}, 



X — Xr- \ 

y-Vr- j 


and 0£*|i<'+ 


1 ( X — Xr+ ^ 

2\K+\ \y-yr+ )' 


Thus lines 4-5 of BDMrightorder is unnecessary for RTq elements. 

For all RTk and BDMk spaces, basis functions are divided into two categories, edge 
based functions and element-based functions. The element based functions are usually 
easy to construct since they are only non-zero in one element, see [1|. So we will only 
discuss edge based basis functions. Both RTk and BDMk spaces have k + 1 basis functions 
on each edge. We only need to functions to span Pk{Eii) on For BDM 2 and i?r 2 , three 
edge basis functions on are 

= A^V-^At, 0^2 = AsAt(V-^Ai - V-^As), and = -XfV-^Xg. 

That is, we use A^, A^, and A^A* to span P 2 on Ei. We can use other choices, for example, 
Legendre polynomials for high order spaces. With this observation, the code developed in 
this paper can be easily adapted to these spaces. 


8.2. Nedelec spaces in 2D. For Nedelec spaces, the 1st and 2nd types Nedelec spaces 
correspond to their H{div) counterparts are RT and BDM spaces, respectively. We also 
only need to discuss the construction of basis functions on edges. By (4.2), the zero- 
moment Ff(curl) edge basis function is 

= A.VAt - AtVA,. 

The linear //(curl) edge basis functions are 

^^f = A,VAi and =-XtV Xg. 


8.3. i/(div) and //(curl) and basis functions in 3D. For a tetrahedral mesh, //(div) 
basis functions are defined on faces. Like the edge structure in this paper, we should 
define a face matrix. If = {zr, Zg, zt}, we need ensure r < s < t, and choose the 
normal direction such that the area of {zr,Zg, zt} is positive. Once this is done, we can do 
similar things as in 2D. The BDMi basis functions in 3D on a face = {zr,Zg,zt} are 

XrVXg X VAt, XgVXt X VXr, and AtVA^ x VA*. 

Its RTq basis function on F is 

At’VAs X VAt Xg\/Xf X Xr “t“ AtVA/. x \/Xg. 

For three dimensionalal //(curl) space, basis functions in barycentric coordinates can be 
found in [12]. 
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